
/****************************************************************************/
/*Copyright (c) 2011, Florent DEVILLE.                                      */
/*All rights reserved.                                                      */
/*                                                                          */
/*Redistribution and use in source and binary forms, with or without        */
/*modification, are permitted provided that the following conditions        */
/*are met:                                                                  */
/*                                                                          */
/* - Redistributions of source code must retain the above copyright         */
/*notice, this list of conditions and the following disclaimer.             */
/* - Redistributions in binary form must reproduce the above                */
/*copyright notice, this list of conditions and the following               */
/*disclaimer in the documentation and/or other materials provided           */
/*with the distribution.                                                    */
/* - The names of its contributors cannot be used to endorse or promote     */
/*products derived from this software without specific prior written        */
/*permission.                                                               */
/* - The source code cannot be used for commercial purposes without         */ 
/*its contributors' permission.                                             */
/*                                                                          */
/*THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS       */
/*"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT         */
/*LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS         */
/*FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE            */
/*COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,       */
/*INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,      */
/*BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;          */
/*LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER          */
/*CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT        */
/*LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN         */
/*ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE           */
/*POSSIBILITY OF SUCH DAMAGE.                                               */
/****************************************************************************/

#include "RTComplex.h"
#include <math.h>

namespace RT
{
	/*Default contructor*/
	RTComplex::RTComplex():m_r(0), m_i(0){}

	/*Constructor*/
	RTComplex::RTComplex(F64 real, F64 ima):m_r(real), m_i(ima)
	{}

	/*Add two complex numbers*/
	RTComplex RTComplex::operator+ (const RTComplex& c)const
	{
		return RTComplex(c.m_r + m_r, c.m_i + m_i);
	}

	/*Substract two complex numbers*/
	RTComplex RTComplex::operator- (const RTComplex& c)const
	{
		return RTComplex(m_r - c.m_r, m_i - c.m_i);
	}

	/*Substract a complex number with a float*/
	RTComplex RTComplex::operator- (F64 f)const
	{
		return RTComplex(m_r - f, m_i);
	}

	/*Unary minus operator*/
	RTComplex RTComplex::operator- ()const
	{
		return RTComplex(-m_r, -m_i);
	}

	/*Multiply two complex numbers*/
	RTComplex RTComplex::operator* (const RTComplex& c)const
	{
		return RTComplex(m_r * c.m_r - m_i * c.m_i, m_r * c.m_i + m_i * c.m_r);
	}

	/*Multiply a complex number with a float*/
	RTComplex RTComplex::operator* (F64 f)const
	{
		return RTComplex(f * m_r, f * m_i);
	}

	/*Division between a complex number and a float*/
	RTComplex RTComplex::operator/ (F64 f)const
	{
		return RTComplex(m_r / f, m_i / f);
	}

	/*Divide a complex number a float and assign the result*/
	void RTComplex::operator/= (F64 f)
	{
		m_r /= f;
		m_i /= f;
	}

	/*Multiply two complex numbers and assign the result*/
	void RTComplex::operator*= (const RTComplex& c)
	{
		m_r = m_r * c.m_r - m_i * c.m_i;
		m_i = m_r * c.m_i + m_i * c.m_r;
	}

	/*Add two complex numbers and assign the result*/
	void RTComplex::operator+= (const RTComplex& c)
	{
		m_r += c.m_r;
		m_i += c.m_i;
	}

	/*Division with two complex numbers*/
	RTComplex RTComplex::operator/ (const RTComplex& c)const
	{
		RTComplex denom = c * c.conjugate();
		F64 fDenom = denom.m_r;

		RTComplex result = *this * c.conjugate();
		result = result / fDenom;
		return result;
	}

	/*Test equality between two complex numbers*/
	bool RTComplex::operator==(const RTComplex& c)const
	{
		if(m_r == c.m_r && m_i == c.m_i)
			return true;

		return false;
	}

	/*Test equality between a complx and a float*/
	bool RTComplex::operator==(F32 f)const
	{
		if(m_r == f && m_i == 0)
			return true;

		return false;
	}

	/*Return the modulus*/
	F64 RTComplex::modulus()const
	{
		return sqrt(m_r * m_r + m_i * m_i);
	}

	/*Return the modulus to the square*/
	F64 RTComplex::modulus2()const 
	{
		return m_r * m_r + m_i * m_i;
	}

	/*Return the square root*/
	RTComplex RTComplex::squareRoot()const
	{
		if(m_i == 0 && m_r == 0)
			return RTComplex();

		if(m_i != 0)
		{
			F64 m = modulus();
			F64 s = m_i / fabs(m_i);
			return RTComplex(sqrt((m + m_r) * 0.5), s * sqrt((m - m_r) * 0.5));
		}
		else if(m_r < 0)
			return RTComplex(0, sqrt(fabs(m_r)));
		else
			return RTComplex(sqrt(m_r), 0);
	}

	/*Return the conjugate*/
	RTComplex RTComplex::conjugate()const
	{
		return RTComplex(m_r, -m_i);
	}

	/*Return the inverse*/
	RTComplex RTComplex::inverse()const
	{
		return RTComplex(m_r / this->modulus2(), -m_i/this->modulus2());
	}

	/*Test if the complex number is a real number*/
	bool RTComplex::isReal()const
	{
		if(m_i == 0)
			return true;

		return false;
	}

	/*Return a stream with the values of the complex number*/
	ostream& operator<<(ostream& s, const RTComplex& c)
	{
		if(c.m_i < 0)
			s << c.m_r << " - i"<<fabs(c.m_i);
		else
			s << c.m_r << " + i"<<c.m_i;
		return s;
	}

}